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Abstract 

In this work , we present the analytic treatment of the 
large order Bessel functions that arise in the Fourier 
Transform (FT) of the Gravitational Wave (GW) signal 
from a pulsar. We outline several strategies which employ 
asymptotic expansions in evaluation of such Bessel func¬ 
tions which also happen to have large argument. Large 
order Bessel functions also arise in the Peters-Mathews 
model of binary inspiralling stars emitting GW and sev¬ 
eral problems in potential scattering theory. Other ap¬ 
plications also arise in a variety of problems in Ap¬ 
plied Mathematics as well as in the Natural Sciences 
and present a challenge for High Performance Comput- 
ing(HPC). 


1. Introduction 

The detection of gravitational waves (GW) from as- 
trophysical sources is one of the most outstanding prob¬ 
lems in experimental gravitation today. Large laser 
interferometric gravitational wave detectors like the 
LIGO, VIRGO, LISA, TAMA 300, GEO 600 and AIGO 
are potentially opening a new window for the study of 
a vast and rich variety of nonlinear curvature phenom¬ 
ena. 

In recent works |T| we have analyzed the Fourier 
transform (FT) of the Doppler shifted GW signal from 
a pulsar with the use of the Plane Wave Expansion 
in Spherical Harmonics (PWESH). Spherical-harmonic 
multipole expansions are used throughout theoretical 
physics. The expansion of a plane wave in spherical har¬ 
monics has a variety of applications not only in quan¬ 
tum mechanics and electromagnetic theory [3] , but also 


in many other areas. A number of researchers have used 
spherical-harmonic expansions for a variety of problems 
in general relativity, including problems where nonlin¬ 
earity shows upgj. The basis states in the PWESH 
expansion form a complete set and facilitate such a 
study. It also turns out that the consequent analysis of 
the Fourier Transform (FT) of the GW signal from a 
pulsar has a very interesting and convenient develop¬ 
ment in terms of the resulting spherical Bessel, gener¬ 
alized hypergeometric function, the Gamma functions 
and the Legendre functions. Both rotational and or¬ 
bital motions of the Earth and spindown of the pul¬ 
sar can be considered in this analysis which happens 
to have a nice analytic representation for the GW sig¬ 
nal in terms of the above special functions. The signal 
can then be studied as a function of a variety of dif¬ 
ferent parameters associated with both the GW pulsar 
signal as well as the orbital and rotational parameters. 
The numerical analysis of this analytical expression for 
the signal offers a challenge for fast and high perfor¬ 
mance parallel computation. The plane wave expansion 
approach was also used by Bruce Allen and Adrian C. 
Ottewill in their study of the correlation of GW sig¬ 
nals from ground-based GW detectors. They use the 
correlation to search for anisotropies from stochastic 
background in terms of the Z,m multipole moments. 
Our PWESH formalism enables a similar study. Re¬ 
cent studies of the Cosmic Microwave Background Ex¬ 
plorer have raised the interesting question of the study 
of very large multipole moments with angular momen¬ 
tum l and its projection to going up to very large val¬ 
ues of l ~ 1000. Such problems warrant an intensive 
analytic study supplemented by numerical and paral¬ 
lel computation. 

Since our FT depends on the Bessel function, a com- 


putational issue arises due to large values of the in¬ 
dex or order n of the function. In the GW form of the 
pulsar, the Doppler shifted orbiting motion gives rise 
to Bessel functions J n {—where 2n fo^ slnt> j s 
large for non-negligible angle 9 as is shown in the fol¬ 
lowing section. Even for sin# ~ the argument is 
large necessitating the consideration of large values of 
n. The motivation of this work, is to extend the anal¬ 
ysis in Watson p] for large index, argument and over¬ 
lapping situations. Meissel f3 has made derivations for 
large order Bessel functions both when the argument 
is smaller than the order and vice versa. The asymp¬ 
totics of these large order Bessel functions are tricky in 
the sense that one runs into so-called “transition” re¬ 
gions where such expansions fail. These regions are val¬ 
ues of the function when the argument is close to the 
given order. As an application, we will address the phe¬ 
nomenological situation of GW signal analysis of large 
order n (which does arise with combinations of l and 
m) and supplement the related computations with the 
presently derived results in a forthcoming paper. 

Captures of stellar-mass compact objects (CO) by 
massive black holes are important capture sources for 
the Laser Interferometer Space Antenna (LISA), the 
space based GW detector due to be launched in about 
a decade [7). Higher Harmonics of the orbital frequency 
of the COs arise in the post Newtonian (PN) capture 
GW model forms and contribute considerably to the 
total signal to noise (S/N) ratio of the waveform. The 
GW form can be decomposed into gravitational multi¬ 
pole moments which are treated in the Fourier analy¬ 
sis of Keplerian eccentric orbits. The radiation depends 
strongly on the orbital eccentricity e, and Bessel func¬ 
tions J„ (ne) are a natural consequence of the analy¬ 
sis. 

The calculation of partial derivatives of the poten¬ 
tial scattering phase shifts which often contain Bessel 
and Legendre functions of large order angular momen¬ 
tum l , with respect to angular momentum arise in a va¬ 
riety of scattering problems in atomic, molecular and 
nuclear physics. In particular, large values of l can arise 
in rainbow, glory and orbit scattering. The analysis in 
our paper should help provide suitable approximations 
for large order and/or argument for the Bessel func¬ 
tions that arise in such problems. 

2. Fourier Transform of the GW signal 

The FT for the GW Doppler shifted pulsar signal 
p.) is given as follows: 
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The angle a is the co-latitude detector angle and an¬ 
gles 9 , (j) are associated with the pulsar source. Here 
c^o 27rfo, uj 0 rb — T or b (T>r& = 365 days, T rE — 1 

day), B orb = 2 (^s. + f , k = ^ IoRe sin(a) 

( R e is the radius of Earth, c is the velocity of light) 
and A = 1.5 x 10 11 meters is the sun-earth distance. 


3. Extensions of Meissel’s and Steepest 
Descent Expansions 

The Bessel function, of the type, J„ (x) obeys the fol¬ 
lowing differential equation 


2 d 2 J v {vz) | dJ v (vz) , 2/1 2 \ j f \ n 
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dz 


where the argument x is parameterized by vz. If a func¬ 
tion u(z) is introduced such that 


J v (yz) = 
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where u(z) is a series in descending powers of v, 
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Substitution of this series and equation ( 8 ) in the dif¬ 
ferential equation (7) yields the following expressions 
for Ui(z ), i = 0...5, 
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Hence, by integrating m, and substituting in Equation 
(8) we arrive at Meissel’s First expansion [5], which is 
valid for the case when the argument is less than the 
order v. We do not list uq, u 7 . u 3 and ug as one can ob¬ 
tain these straightforwardly from their respective inte¬ 
grals shown below. These results are expressed as, 

T t ^ _ {vz) v exp(i/Vl - z 2 ) exp(-K) , 

A ’ e v T{v + 1)(1 — z 2 ) 4 / 4 [l + vT - ”? 2 ]'' 

where, 

V v = V 1 + V 2 + V 3 + V 4 + V 5 + V 6 + V 7 + V 8 + ... (11) 
and, 
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Hence we have actually increased Meissel’s analysis by 
two orders. Using symbolic packages these orders were 
computed and higher terms should pose no problem if 
the application requires higher accuracy. 

For the case when the argument is larger than the 
index, Meissel used the parametrization z = sec (3 0, 
and we shall term it as his Second expansion. Hence, 
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where P„ is given as, 
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and, 
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It should be remarked that we disagree with Meissel’s 
result for P 3 in the last four terms. However, we ob¬ 
tain perfect agreement with the rest of his results |Bj. 

We have improved on his result by using V?, V$ to ob¬ 
tain Pa and Q 4 . Hence, we have increased the accuracy 
of this expansion by at least one order from Meissel’s 
earlier result. Again, higher order results are easily ob¬ 
tainable and are available if needed. 

In Figures 1 and 2 we have plotted these expan¬ 
sions in the regions they are expected to fail. These 
are the so called “transition” regions, where each ex¬ 
pansion approaches a singularity (as the order equals 
the argument). For the computationally motivated (we 
can compute exact values of Bessel functions with ease) 
case of the v = 300, we note the following. Fig.l indi¬ 
cates the onset of breakdown in the First expansion 
for argument values around and larger than 290. Sim¬ 
ilarly, Figure 2, indicates a similar breakdown starting 
around the values 300 and persisting till 310. Hence, 
the values outside these regions of breakdown or tran¬ 
sition regions are well covered by Meissel’s expansions. 

However, the issue as to deal with these regions need to 
be addressed via separate methods, which will be ad¬ 
dressed in more detail in Section IV. The CPU time 
for these approximations was less than 0.01 seconds 
per value on a 2.4 GHz Pentium IV processor running 
MAPLE version 9. The “exact” MAPLE solver took 
somewhere between 0.03 to 0.08 seconds to compute 
each value. Clearly, there is a lot more computational 
speed in using a few terms present in these expansions. 

As an application, it should be noted that values of 
this order are applicable to the Peters-Mathews model 
of gravitational radiation from binary inspiralling stars 
0 


function graphed for argument x and order v = 300 
near the transition region. Solid line indicates actual 
Bessel function values and circles indicate values given 
by the expansion. 



Figure 2: Meissel’s Second expansion and actual Bessel 
function graphed for argument x and order v = 300 
near the transition region. Solid line indicates actual 
Bessel function values and circles indicate values given 
by the expansion. 





















For the case when the argument equals the index, 
we extend Meissel’s Third expansion 0 by two orders 
as follows: 
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where the terms, A m (m = 0,1, 2, ..7), are given by, 
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We observe that inclusion of the higher order terms 
leads to 10 decimal accuracy compared to actual val¬ 
ues of large order Bessel functions. 

The method of steepest descents was employed by 
Debye in For the case when the argument is less 
than the order, he obtained, 
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Following this method, we have computed two higher 
orders A 3 and A 4 , using symbolic computation. 

For the case when the argument is larger than the 
order, Debye obtains the following expansion: 
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Again, we have extended Debye’s result by two 
higher orders by obtaining A 3 and A 4 . However, due to 
the nature of this method we could not obtain reliable 
results that spanned in a generally predictable direc¬ 
tion. Accuracy was limited to the region of the station¬ 
ary phase as expected and hence, we recommend Meis- 
sel’s expansions to be more reliable (except of course 
in the “transition” region) than the method of steep¬ 
est descent. 

4. Transitional regions: Contour Inte¬ 
gration and extension of e expansion 

To address the issues related to computation for 
large order Bessel functions in the transition regions 
we present two methods that are geared to work in 
such domains. Firstly, we present the results by Wat¬ 
son, 0. For the case of the argument being less than 
the order, he obtained via use of contour integration, 
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where 9 X < 1. Similarly, for the case when the argu¬ 
ment is greater than the order, he derived the follow¬ 
ing: 
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where 0 2 < 1 and the argument for the Bessel func¬ 
tions J ± i is )r tan 3 / 3. The great advantage of these for¬ 
mulae is that they have error bounds given. However, 
these extensions are not trivial as this involves solv¬ 
ing extensions to Airy-type integrals, for which we do 
not presently have closed form answers. The other is¬ 
sue with these formulae is that they are themselves 
given in fractional Bessel function form which would 
pose computational problems once the arguments in¬ 
volved are large. 

On the other hand, Debye [H], introduced, what we 
will term as “e expansion”. The idea is motivated by 
introducing a small parameter e, such that v = z( 1 — 
e), where v denotes the order and z is the argument of 
the Bessel function. 
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We have extended this analysis by 5 orders and the 
terms B m (ez), m = 0,1, 2, ..15, are given as, 
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Terms B 3m _i , to = 1,2... do not contribute in eqn. 
(21) due to the periodicity of the sine function. With 
symbolic computation, we are able to generate higher 
orders if needed. 

To illustrate the applicability and issues of both 
these methods to the transition region, we present Fig¬ 
ures 3 and 4, which are plotted for the problematic re¬ 
gions (when the order is v = 300) in Figures 1 and 2. 
Both methods show remarkable ability in capturing the 
functions in the domains of interest. In Figure 3, the 
e expansion starts working at values at 286 and Wat¬ 
son’s formula works to even a larger domain. Similarly, 
in Figure 4, both the methods indicate success in re¬ 
gions where Meissel’s expansions fail. This starts at val¬ 
ues of the argument, and works up to x = 316 for the e 
expansion whereas, again, the domain of Watson’s for¬ 
mula is much greater. The reasons for lesser range of 
the e expansion can be attributed to the fact that it 
is a power series compared to Watson’s formula which 
actually depends on fractional Bessel functions them¬ 
selves. Further, the e expansion depends crucially on 
the size of the parameter, which is connected with the 
order one is working with. However, the reason why we 
will persist with this method is that it will be more ap¬ 
plicable when the argument of the Bessel function is 
quite large. 

To illustrate the type of values a GW pulsar FT 
would require, we present Figures 5 and 6. Here, we 
choose a very large order (yet realistic phenomeno¬ 
logically) for the Bessel function, which is 1 million. 












































Figure 3: Comparison of e expansion and Watson’s for¬ 
mulae for argument x < 300 and order v = 300.Solid 
line indicates exact Bessel function values, diamonds 
represent e expansion and circles indicate values given 
by Watson’s formula. 



Figure 4: Comparison of e expansion and Watson’s for¬ 
mula in the transition region for argument x > 300 and 
order v = 300. Solid line indicates exact Bessel func¬ 
tion values, diamonds represent e expansion and circles 
indicate values given by Watson’s formula. 


Also, in such a scenario, we would be looking at values 
greater than one million, hence Meissel’s second expan¬ 
sion along with the appropriate Watson’s formula (eq. 
20) will be put to use. We were not able to make exact 
comparison, obviously due to massive computer times 
required. In this regard, the problem of “exact” Bessel 
functions presents a genuine challenge to SHARCNET 
(Shared Hierarchical Academic Research Cluster Net¬ 
work) and HPC in general. In Figure 5, we observe 
strong evidence that the proposed asymptotic expan¬ 
sions are appropriate for GW signal analysis. Here, we 
note the transition region starting at values of the ar¬ 
gument at 1,000,000 and going up to 1,000,200. In this 
region, both the e expansion and Watson’s formula al¬ 
most coincide with each other. As usual, the e expan¬ 
sion breaks down earlier, however, all three methods co¬ 
incide in a certain region indicating that we have con¬ 
sistent methods that work for values relevant to GW 
analysis. Meissel’s expansion is fairly easy to imple¬ 
ment computationally and indicates good stability for 
rather large values of the argument. This is illustrated 
in Figure 6, where we plot this expansion for values 
ranging from 1,000,200 to 32,500,000, which are rele¬ 
vant for GW phenomenology. This appears as a black 
band and is a continuous function which indicates oscil¬ 
lations tightly bunched together. It is noteworthy that 
the method is stable and shows consistent behaviour 
over an extreme range of values for the argument. The 
CPU time consumed by each of the points, on the aver¬ 
age took less than 0.01 seconds on MAPLE. The Bessel 
utility in MAPLE crashed repeatedly after 15-30 min¬ 
utes on the same system described above. It should be 
remarked that Watson’s formula lacks in this capac¬ 
ity as it depends on fractional Bessel functions itself, 
which will provide computational challenge for such 
values. A detailed analysis regarding computational ad¬ 
vantage over exact computation will be addressed in a 
later work. It is aimed to not only address the ques¬ 
tion of GW analysis but will deal with general compu¬ 
tational issues regarding large order Bessel functions. 


5. Conclusion 

In this present work, we have given an extended 
asymptotic analysis for the large order and argument 
Bessel functions. This analytically improves the earlier 
pioneering works of Meissel, Airey, Debye and Watson. 
These extensions should be of possible use not only 
in GW signal analysis, but also in a variety of prob¬ 
lems in Engineering and the Sciences where the ubiq¬ 
uitous Bessel functions are encountered. 
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